# 2023.9.15
library(Seurat)
library(tidyseurat)
library(tidySingleCellExperiment)
library(ggpubr)
library(DropletUtils)
source("00_util_scripts/mod_seurat.R")

# load seurat objects ----------
sobj_sg <- read_rds('CRC-I/data/seekgene/crc-starsolo-annotated.rds')

sobj_10x <- read_rds('CRC-I/data/zy2020_tumor10x.rds')

sobj_guo <- read_rds("CRC-I/data/guo2021sided/guo2021_immune.rds")

sobj_liv <- read_rds('CRC-I/data/li_2021_liver/li2021crc.rds')

# merge 4 datasets ----------
sobj_guo$dataset <- 'guo2021'
sobj_liv$dataset <- 'li2021'
sobj_sg$dataset <- 'liu2023'
sobj_10x$dataset <- 'zhang2020'

gene_merged <- tibble(gene = rownames(sobj_10x)) |>
  filter(gene %in% rownames(sobj_sg),
         gene %in% rownames(sobj_liv),
         gene %in% rownames(sobj_guo)) |>
  pull(gene)

gene_merged3 <- tibble(gene = rownames(sobj_sg)) |>
  filter(gene %in% rownames(sobj_liv),
         gene %in% rownames(sobj_guo)) |>
  pull(gene)

# this can take ~1 min
sobj_c4 <- c(sobj_sg, sobj_10x, sobj_liv, sobj_guo) |>
  reduce(merge) |>
  subset(features = gene_merged)

sobj_c4 <- sobj_c4 |>
  select(.cell, orig.ident:nFeature_RNA,
         genotype, dataset, hpca_main)

sobj_c4 <- sobj_c4 |>
  quick_process_seurat(batch = c('orig.ident','dataset'),
                       pcs = 25,
                       skip_norm = TRUE)

sobj_c4 |> DimPlot(cols = DiscretePalette(36))
sobj_c4 |> DimPlot(split.by = 'dataset',ncol = 2)

## manually annotate ----------
sobj_c4 |> DotPlot(features = seurat_markers) + RotatedAxis()

### marker for ambiguous clusters -------
ambi_markers <- c(11,16:21) |>
  map(\(x)FindMarkers(sobj_c4, ident.1 = x, only.pos = TRUE, logfc.threshold = 2)) |>
  map(rownames_to_column, 'gene') |>
  set_names(c(11,16:21)) |>
  list_rbind(names_to = 'cluster')

ambi_markers |>
  filter(cluster == 16) |>
  slice_max(avg_log2FC, n = 10)

no_other_b_11 <- sobj_c4 |>
  filter(seurat_clusters != 3 & seurat_clusters != 4) |>
  FindMarkers(ident.1 = 11, only.pos = TRUE, logfc.threshold = 1.5)

# comparing dotplot and dimplot to manual annotate merge4, do not delete this!
sobj_c4 <- sobj_c4 |> mutate(manual_main = case_match(
  as.numeric(seurat_clusters),
  c(3,4) ~ 'B',
  c(1,5,10,12) ~ 'CD4T',
  c(2,6,18) ~ 'CD8T',
  c(7,9,14:17) ~ 'Myeloid',
  8 ~ 'Mast',
  13 ~ 'NK',
  11 ~ 'Dividing lymphocyte',
  c(19,21) ~ 'Epithelia'))

sobj_c4 <- sobj_c4 |>
  filter(!is.na(manual_main))

sobj_c4 |> DimPlot(group.by = 'manual_main', label = TRUE,
                   label.box = TRUE) +
  ggtitle('Major immune cells')

sobj_c4 |> DimPlot(group.by = 'manual_main', split.by = 'dataset',ncol = 2) +
  ggtitle('Major immune cells')

## auto annotate ------
sobj_c4 |> DimPlot(group.by = 'hpca_main', label = TRUE,
                   label.box = TRUE)

monaco <- celldex::MonacoImmuneData()

sobj_c4 |>
  filter(seurat_clusters %in% c(3,4,11)) |>
  mark_cell_type_singler(monaco, fine_label = TRUE, new_label = 'monaco_fine')

sobj_c4 <- sobj_c4 |>
  mark_cell_type_singler(hpca, new_label = 'hpca_main_c4')

sobj_c3 <- sobj_c3 |>
  mark_cell_type_singler(hpca, new_label = 'hpca_main_c3')

sobj_c4 |> DimPlot(group.by = 'hpca_main_c4',
                   cols = DiscretePalette(10))+
  ggtitle('Major immune cells')

sobj_c3 |> DimPlot(group.by = 'hpca_main_c3',
                   cols = DiscretePalette(10))+
  ggtitle('Major immune cells')

sobj_c4 |> DimPlot(group.by = 'hpca_main_c4', split.by = 'dataset',ncol = 2,
                   cols = DiscretePalette(10)) +
  ggtitle('Major immune cells')

### eliminate ambiguous & non-immune cells -------
marker11 <- FindMarkers(sobj_c4, ident.1 = 11, only.pos = TRUE, logfc.threshold = 1)

marker11 |>
  as_tibble(rownames = 'gene') |>
  arrange(desc(avg_log2FC))

DotPlot(sobj_c3,features = cc.genes.updated.2019)

sobj_c4 |> write_rds('CRC-I/data/crc_merge4_immune.rds')

sobj_c4 |> mutate(genotype = ifelse(str_detect(orig.ident, 'right'), 'II', genotype)) -> sobj_c4

sobj_c4 <- read_rds('CRC-I/data/crc_merge4_immune.rds')

sobj_c4@meta.data |>
  as_tibble(rownames = '.cell') |>
  write_csv('CRC-I/results/crc_merge4_meta.csv')

mast_tt_deg <- sobj_c4 |>
  filter(manual_main == 'Mast') |>
  FindMarkers(ident.1 = 'TT', group.by = 'genotype')

nk_tt_deg <- sobj_c4 |>
  filter(manual_main == 'NK') |>
  FindMarkers(ident.1 = 'TT', group.by = 'genotype')

main_tt_deg <-sobj_c4$manual_main |>
  unique() |>
  map(\(x)sobj_c4 |> filter(manual_main == x) |>
        FindMarkers(ident.1 = 'TT', group.by = 'genotype')) |>
  map(rownames_to_column, 'gene') |>
  set_names(unique(sobj_c4$manual_main)) |>
  list_rbind(names_to = 'cluster')

main_tt_deg |> write_csv('CRC-I/results/merge4_main_tt_deg.csv')

# subcluster ------
sobj.cd8 <- sobj |>
  filter(str_detect(zhang2020_main, 'CD8'))

sobj.cd8 <- sobj.cd8 |>
  GetAssay() |>
  CreateSeuratObject(meta.data = sobj.cd8@meta.data)

sobj.cd8 <- sobj.cd8 |> quick_process_seurat()

sobj.cd8 <- sobj.cd8 |>
  mark_cell_type_singler(
    ref = filter(sce_10x, str_detect(label.main,'CD8')),
    sc_ref = TRUE,
    fine_label = TRUE,
    new_label = "zhang2020_fine"
  )

## examine FCGR2B expression ---------
fcgr2b_expr <- sobj |>
  get_abundance_sc_long('FCGR2B')

fcgr2b_expr <- sobj |> as_tibble() |>
  select(.cell, zhang2020_main) |>
  left_join(fcgr2b_expr)

fcgr2b_expr |>
  ggplot(aes(zhang2020_main, .abundance_RNA, color = zhang2020_main)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white', linewidth = 2) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_cl_boot', width = .3) +
  ggpubr::theme_pubr() +
  labs(x = 'Cell type', y = 'Normalized expression', color = 'Cell type') +
  NoLegend()

g1 <- last_plot()

fcgr2b_expr |>
  group_by(zhang2020_main) |>
  mutate(wt = 1/n()) |>
  filter(.abundance_RNA > 0) |>
  count(zhang2020_main, wt = wt) |>
  ggplot(aes(zhang2020_main, n, fill = zhang2020_main)) +
  geom_col() +
  ggpubr::theme_pubr() +
  labs(x = 'Cell type', y = 'Fraction of FCGR2B+ cells', color = 'Cell type') +
  NoLegend()

g2 <- last_plot()

g1 + g2 + patchwork::plot_annotation(title = 'CRC TME immune cell expression of FCGR2B (Guo, 2021, JCI Insight)')

### only cd8 -----
fcgr2b_expr <- sobj.cd8 |>
  get_abundance_sc_long('FCGR2B')

fcgr2b_expr <- sobj.cd8 |> as_tibble() |>
  select(.cell, zhang2020_main) |>
  left_join(fcgr2b_expr)

fcgr2b_expr |>
  ggplot(aes(zhang2020_main, .abundance_RNA, color = zhang2020_main)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white', linewidth = 2) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_cl_boot', width = .3) +
  ggpubr::theme_pubr() +
  labs(x = 'Cell type', y = 'Normalized expression', color = 'Cell type') +
  NoLegend() +
  coord_flip() +
  scale_color_brewer(type = 'qual')

g1 <- last_plot()

fcgr2b_expr |>
  group_by(zhang2020_main) |>
  mutate(wt = 1/n()) |>
  filter(.abundance_RNA > 0) |>
  count(zhang2020_main, wt = wt) |>
  ungroup() |>
  add_case(zhang2020_main = 'LAYN-CD8-Tex', n = 0) |>
  ggplot(aes(zhang2020_main, n, fill = zhang2020_main)) +
  geom_col() +
  ggpubr::theme_pubr() +
  labs(x = 'Cell type', y = 'Fraction of FCGR2B+ cells', color = 'Cell type') +
  NoLegend()+
  coord_flip() +
  scale_fill_brewer(type = 'qual')

g2 <- last_plot()

g1 + g2 + patchwork::plot_annotation(title = 'CRC TME CD8 T cell expression of FCGR2B (Guo, 2021, JCI Insight)')

sobj |> write_rds("CRC-I/data/guo2021sided/guo2021_ident.rds")

